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ABSTRACT 

In hierarchical clustering, galaxy clusters accrete mass through the aggregation of 
smaller systems. Thus, the velocity field of the infall regions of clusters contains signif¬ 
icant random motion superimposed on radial infall. Because the purely spherical infall 
model does not predict the amplitude of the velocity field correctly, methods estimat¬ 
ing the cosmological density parameter Hq based on this model yield unreliable biased 
results. In fact, the amplitude of the velocity field depends on local dynamics and only 
very weakly on the global properties of the universe. 

We use A-body simulations of flat and open universes to show that the amplitude 
of the velocity field of the infall regions of dark matter halos is a direct measure of the 
escape velocity within these regions. We can use this amplitude to estimate the mass 
of dark matter halos within a few megaparsecs from the halo center. In this region 
dynamical equilibrium assumptions do not hold. The method yields a mass estimate 
with better than 30% accuracy. If galaxies trace the velocity field of the infall regions 
of clusters reliably, this method provides a straightforward way to estimate the amount 
of mass surrounding rich galaxy clusters from redshift data alone. 


Subject headings: dark matter — galaxies: clusters: general — gravitation — methods: 
numerical 


1. INTRODUCTION 


The linear theory of density perturbations shows that a spherically symmetric mass concen¬ 
tration in an expanding universe induces a radial peculiar velocity held in the surrounding region 

Wec(^) 

where 6{r) is the average spherical mass overdensity within the radius r, and Hq and Hq are the 
Hubble constant^ and the cosmological density parameter at the present time, respectively. We 


^Present address: Max-Planck-Institut fiir Astrophysik, Karl-Schwarzschild-Str. 1, 85740, Garching bei Miinchen, 
Germany 

^We use h — 0.5 throughout, where Hq = lOOh km s“^ Mpc“^. 
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can estimate the galaxy number overdensity 6g = b5 where b is the bias parameter. Thus, if we can 
measure Upec we can estimate (3 = QQ^/b. Several authors have applied this method to the Local 


Supercluster and to the Virgo cluster with mixed results (see e.g. the review by Davis fc Peebles 
see also ptrauss &: Willick 1995). Along with the uncertainties in the determination of the 


1983 


galaxy number overdensity and the peculiar velocity themselves, external tidal shear can strongly 
affect the velocity held (e.g. Hoffman 1986 ; Eisenstein Sz Loeb 1995 ; Bond Sz Myers 1996| ); thus we 
cannot obtain reliable results unless we sample the velocity held within the full three-dimensional 
region around the cluster (Villumsen Sz Davis 198^). 


In redshift space, spherical infall conhnes galaxies around clusters within caustics, surfaces 
with a characteristic “trumpet” shape ( Kaiser 19871) . For these non-linear regions, we must replace 
the linear regime equation (||) with the exact solution of the equation of motion of shells in an 
overdense region within the expanding universe (e.g. Silk 197^. However, the dependence of the 


peculiar velocity on (3 and 5g is still approximately separable ( Regos &: Geller 1989 ). Thus, provided 
we can determine the location of the caustics, we can still estimate (3 through the measurement of 
6g for a rich galaxy cluster ( Regos &: Geller 1989 ). Recently, [Regos (1996)| suggested the application 
of gravitational leasing (e.g. Tyson, Valdes, k, Wenk 1990|; [Kaiser &: Squires 1993|; [Kaiser, Squires 


[&: Broadhurst 1995[ ; Bonnet &: Mellier 1995[ ) to determine the mass overdensity <5 and therefore the 
unbiased value of Dq- 


Standard inflationary cosmologies (e.g. Peacock 1996[) predict that the primordial density field 
is a Gaussian random field. These initial conditions can lead to either top-down or bottom-up 
scenarios for the formation of cosmic structures. The hierarchical clustering in the bottom-up 
scenarios, with large systems forming by aggregation of smaller ones, currently represents the most 
successful framework of structure formation theories as it seems to be able to reproduce many, 
although not all, properties of the real universe. The top-down scenarios are far less succesfull (see 
e.g. [Padmanabhan 1993 for a general discussion). 


The general formation process of clusters in hierarchical scenarios differs substantially from 
the one described by the idealized spherical model; here clusters form through the infall of smooth 
spherical shells onto an initial density peak ( Punn &: Gott 1972[ ). Despite its idealized nature, 
the spherical model appears to represent a reasonable description of the collapse of high peaks in 


a random Gaussian density field ( Bernardeau 1994 ), and predicts density and velocity profiles of 
hnal systems in reasonably good agreement with V-body simulations of dark matter halo formation 


when the power spectrum has an effective spectral index n > — 1 ( Zaroubi, Naim, &: Hoffman 1996 ). 
These authors suggest that particle ranks in binding energy are conserved during the formation 
process. This conservation is responsible for the agreement (see also Hoffman 1988[ ; puinn Sz Zurek 
1988; [Zaroubi Sz Hoffman 1991^ ). In fact, for spectral indices n < —1 the formation process 
violent (Lynden-Bell 19^), energy ranks are not conserved, and the agreement breaks down. 


IS more 


Thus, the agreement in energy space between the spherical infall model and hierarchical clus¬ 
tering for a limited range of spectral index, n, allows correct predictions about the final state of 
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the halo. However, the spherical infall model does not describe the evolution of the outer regions 
of systems in conhguration space. A/^-body simulations of flat universes ( 

Haarlem fc van de Weygaert 199^ ) show that the velocity field of cluster surroundings is very differ¬ 
ent from the predictions of spherical infall because (1) recent mergers increase the particle kinetic 
energy, and (2) the presence of substructures makes the velocity profile irregular. In other words, 
random motions obscure the infall information. If the same problem arises in real galaxy clusters, 
estimates of (5 (or flo) based on spherical infall are systematic overestimates. 

Here, we suggest a unifying explanation for the amplitude of the caustics. We use Wbody 
simulations of flat and open universes to show that the escape velocity around dark matter halos 
determines the amplitude of the caustics. Thus, the local dynamics of the halo, not the global 
properties of the universe, dominates the amplitude of the velocity held around clusters. 

We also show that under reasonable hypotheses about the density held outside the virialized 
region, the velocity held around a cluster provides an estimate of the mass enclosed within a few 
megaparsecs of the cluster center. This method is particularly interesting because, in this region, 
the equilibrium assumptions underlying usual mass estimation methods do not hold. 

In Sect. 2 we review the main equations of the spherical infall model. In Sect. 3 we derive the 
alternative expression for the amplitude of the caustics and we outline the mass estimation method. 
In Sect. 4 we compare both the expression for the amplitude of the caustics and the spherical infall 
model with the results of Wbody simulations. We also use these Wbody simulations to test our 
new mass estimation method. 


van Haarlem 1992; 


van 


2. THE SPHERICAL INFALL MODEL 


We now review the main points of the infall model for a spherical perturbation in the case of 
open and flat universes. We then define the infall region and the caustic surfaces surrounding the 
perturbation. 

Consider a spherical perturbation in an expanding universe described by the average overden¬ 
sity profile 5i{< r*) within the physical radius rj at time a* <C 1, where a is the cosmological scale 
factor, a = 1 at the present time, and <?; 1 for any radius r*. The evolution of each shell in the 
absence of shell crossings is described by the equation of motion dFr/dt^ = —GM/r'^ + Ar/3, where 
M is the mass within r, A is the cosmological constant, and G is the gravitational constant. The 
first integral of the equation of motion is 

2 \dt J r 6 

where E is the total energy of the shell. In order to consider the collapse independent of the 
primordial density profile, we can think in terms of the local scale factor x = raijri. Note that at 
time tti, X = Qi for all the shells. However, at later times, if we assume that 6i is a monotonically 


( 2 ) 
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decreasing function of r*, the degeneracy between x and r* disappears, and we have a biunique 
correspondence between x and With this change of variable, equation (|^ becomes 

x^ — — — = e (3) 

X 

where the dot indicates a derivative with respect to r = ^ao = and 

a{6i) = no(l + 5i) (4) 

= 1 “ ^Ao ~ ^0 • (^) 


In deriving equation (^) for the effective energy e = {2E/HQ){ai/ri)‘^, we set the initial condi¬ 
tions at time ai ~ 10“^, when linear theory in a matter-dominated universe is a good description 
of the growth of perturbations. In other words, we assume that at this epoch the decaying modes 
of the perturbation have decayed away and the peculiar velocity field Vi has reached the amplitude 
predicted by linear theory, Vi/Hin = where Hi is the Hubble constant and Hoi is the 


cosmological density parameter at time a* (see e.g. Lilje &: Lahav 1991). Thus, to first order in 
5i, the peculiar velocity field contributes the term —2Q.Q^^f5i/‘iai to the effective total energy. 
In equation (^) we have further assumed Qoi = + (1 — ~ ^Ao)0‘i + — 1) when 

Qi ~ 10“^ <C 1. 


e < 


( 6 ) 


When Hao < 0, the shell always collapses. Otherwise, for collapse to occur, the effective energy 
e must satisfy 

-30Ao(a/2HAo)"/^ ^^Ao > 0 
0 Hao = 0. 

The smallest positive root (or the only positive root when Oao < 0) of the cubic equation OAoa^^ + 
ex + a = ^ gives the maximum expansion scale factor Xmax(<5i) of the shell with internal overdensity 
5i. Thus, we may write the collapse time 


’^coll(^*) — 2 


dy 


lo x{y,6i)' 

Finally, we can write the equation describing the evolution of x{5i,T) implicitly 


r = 


I dy/x{y, 6i) 0 < r < rcoii/2 

1 'rcoii((5i) - dy/x{y, Si) rcoii/2 < r < TcoH- 


(7) 


( 8 ) 


In terms of the local scale factor x, we next determine the caustics which depend only on the 
cosmological parameters and not on the overdensity profile. Consider the spherical perturbation 
within space X = (x,9,ip), where 9 and (p are the azimuthal and longitudinal angle, respectively. 
In the following, we use the terms “distance” and “velocity” to indicate the quantities x and x, 
respectively. Recall that x = vaijri. Thus, we need the initial average overdensity profile Si{< r*) 
to transform space X into configuration space R = {r,9,ip), and to obtain the physical distance r 
and the physical velocity v = dr/dt. 
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Now, suppose we project the spherical perturbation within space X onto the plane 6 = 7 r/ 2 . 
We can measure the projected distance x± from the origin, the longitudinal angle ip and the velocity 
component X 3 along the axis 0 = 0. The relations between the new and the old variables are 

x_L = x((5i, r) sin 0, 0:3 = x(x, <5^) cos 0, ip = ip. (9) 


Suppose we observe the spherical perturbation at a time, r. From the equations Tcoii((5f) = r 
and rcoii(5™) = 2 r (eq. [^) we can determine the primordial overdensity of the shell which 
has just collapsed and the (5™ of the shell which is at maximum expansion at time r. Thus, the 
local scale factor of the perturbation at turnaround is xta = To determine the virial 

scale factor at time r we need to consider the virial theorem (see e.g. [Lahav et al. 19M ). We 
have e = —a/ 2 x^ — 2 Q\qxI in virial equilibrium, and e = —ajxm — same shell at 

maximum expansion. Energy conservation implies 


2riil)^ — (2 + r])^l; + 1 = 0 


( 10 ) 


where 'll: = x^^lxm and 7 / = 2Vt.\Qx‘^/a. Note that we can also write 7 / = A/dvrCpm, where pm is the 
mean density of the perturbation within = ViXm/o-i- The condition for collapse is 7 / < 1, and 
we are interested in the solution < 1 of equation ®) with r\ a free parameter. We have 


{ 2pi/3cos[arccos((7/p)/3] V < Po 

(g + Ai/2)V3_g(g_ AV2)1/3 rjo<v<0 

1/2 7/ = 0 

2pi/3 cos[arccos(( 7 /p )/3 + 47 r/ 3 ] 0 < 7 ? < 1 


( 11 ) 


where p = [(2 + r])/ 6 r]]^^'^, q = — 1 / 47 /, A = and t/q — —6.427 is the only real solution of 

the cubic equation A = 0 . 


Now we can write the virial scale factor Xvir at time r, Xvir = where 7/ = 

2nAo37max('^i^)/llo(l + ^i)- We define the infall region of the perturbation at time r as the re¬ 
gion within the scale factor range (xvir,a:ta)- Note that Xvir and xta are the virial scale factor and 
the turnaround scale factor of two different shells of the perturbation. 


At hxed x± G (xvir, a^ta), the velocity component X 3 depends only on 0, ±3 = x{x^/ sin 0,5^(0)] cos 0 
(see eqs. and |^), where 5^(0) satisfies the implicit equation (see eq. |p) 


T = rcoii[5i(0)] - [ 
Jo 


xj_/ sin 6 


dy 


x[v,di{e)\ 


( 12 ) 


Thus, the extrema of ±3 occur at 0 = 0max(a;±, t), and 0 = vr — 0 max (a; 1 ,r) which satisfy the 
equation dx 3 /d 9 = 0. Note that 0max is almost independent of the cosmology. 

In the plane (xj_,X 3 ) the perturbation fills only the region within the curves 


a' 3 (x_L,T) — X {x[<5j (0iiiax) ) x]} cos 0max 


(13) 
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which determine the caustics. For any x±, the maximum difference in the velocity component X 3 
is the difference between the velocities on the caustics, Ax 3 (x^,r) = 2 |i; 3 (x_L, t)|. 

To switch from space X to configuration space i?, we must take the initial average overdensity 
profile 5i{< ri) into account. However, at any time, the physical radius is r = xvijai and the radial 
velocity is v/Hor = x/x. The peculiar velocity is v^ec/H^r = v/Hqt — H/Hq, where H is the 
Hubble constant at time r. Thus, for the velocity difference on the caustics, we have Ax^/x^ = 
Avs/H qvx, where r± is the projected physical radius along the line of sight. Therefore, if we 
define the amplitude of the velocity field = AU 3 / 2 , we can see that the observable quantity 
Avia^rx coincides with the quantity Ax 3 / 2 x_l we have in space X. Fig. 1 shows the profile of 
the amplitude of the velocity field vs. the local scale factor x± of the projected physical radius r± 
for different values of [Ho,Hao] at the present time. We have r^/rvir = {x^/x^\-^){r^/r^\y:)a=ai and 
a monotonically decreasing 5i{< rj); thus we must stretch these curves out according to 5i{< rj) in 
order to obtain the observable profiles. Note that xta/a^vir ^ 2^/^ is almost constant for Hao = 0 
universes. In particular, we have rta/^vir = (a^ta/a^vir)(Ha/?’vir)a=ai ^ 3 in the absence of shell 
crossings. 

Fig. 1 shows that the caustic amplitude depends on [Hq, Hao]) as expected. However, for fixed 
x±_ the dependence is weak: (91nx/51nHo ~ 0.3. We recover the usual 0.6 logarithmic dependence, 
when we switch to physical distances r±. 


3. ESCAPE VELOCITY AND MASS ESTIMATE 


In redshift space, the infall regions of halos in A-body experiments resemble the characteristic 
trumpet shapes expected in the spherical infall model (e.g. van Haarlem &: van de Weygaert 1993 ; 
ling &: Horner 1996| ). However, the amplitude of these caustics, namely the difference between the 
maximum and minimum line-of-sight velocity at projected distance r±_ from the halo center, usually 
exceeds the amplitude predicted by the model. The larger amplitude occurs because of random 
motions in the infall region. 


In §3.1, we derive an expression for the amplitude of the velocity field which takes random 
motions into account. This expression indeed replaces the prediction of the spherical infall model. 
Moreover, this expression suggests a method of estimating the mass enclosed within a distance r± 
from the halo center, where r± extends well beyond the region of virial equilibrium. We derive this 
expression in §3.2. In §4, we compare both expressions with A-body simulations of dark matter 
halos. 






3.1. Amplitude of the velocity field 


Consider a dark halo in its center of mass reference frame. Consider the velocity anisotropy 
parameter at position r from the halo center 


/3(r) = 1 - 



(14) 


where Vr is the radial component of the particle velocity v, = v'f + v'^, and (•) indicates an 
average over the velocities of all particles within the volume d^r centered on position r. Within 
the virialized region r < rg, where is some virial radius (see §4.2.1), particle motion is mainly 
random and /9 ~ 0. For r > r^, orbits become more radial and (3 increases up to ~ 0.8 (see §4.2.1). 
Within these regions, random motions still contribute significantly to the velocity field, because of 
the intrinsic stochastic accretion process in hierarchical clustering. The non-zero (n^) at large radii 
is responsible for the general increase of the amplitude of the caustics relative to pure spherical 
infall. 


In order to write down an expression for the amplitude of the velocity field, we note that 
the presence of random motions in the outskirts of the halo implies that encounters can give a 
particle enough energy to escape the gravitational field of the halo. Thus, if the random motion is 
significant, we expect that the amplitude of the velocity field at a given distance r is the line-of-sight 
component of the escape velocity, not the smaller velocity induced solely by spherically symmetric 
infall. 

The escape velocity is Ugscl'^) = —2i^(r): the caustic amplitude is thus a measure of the 
gravitational potential For the sake of simplicity, consider a spherical system. We have 

^esc(^) = “2(/)(r) = ''") _|_ f p(^x)xdx (15) 

T Jr 

where M(< r) is the total mass within r and p{r) is the system density profile. Equation ( |l5|) holds 
regardless of the stability of the system. 

In equation (|I5| ) v^sc is the full three-dimensional escape velocity. Observations of real clusters 
will provide only the component ui.o.s. of the escape velocity along the line-of-sight at projected 
distances r_\_- The velocity field within the infall region is not isotropic, as indicated by the large 
value of the anisotropy parameter, (3 ^ 0.5. If the tangential component vt[r±) of the velocity field 
is isotropic, we have (uf) = ),2 or 

(^L(^’±)) = (^L.s.(^±)) ^^ 

In §4.2.2, we will see that equation (|l^ indeed describes the amplitude of the velocity field of dark 
matter halos out to a few times the virial radius, r^. 


^We are assuming that Uesc(r) is a non-increasing function of r. This behavior is valid for density profiles steeper 






3.2. Mass estimate 


We now argue that the measure of the escape velocity along the line of sight, i.e. the measure of 
the amplitude of the caustics, may provide a method to estimate M(< r). In principle, knowledge 
of Vesc readily yields an estimate of the mass within r from equation (^), 

fjip' 

2GM(< r) = (17) 

dr 

However, is the product of two functions (eq. |^). The anisotropy parameter (5{r±) is 
unknown, but suppose we can model it. The only measurable function ^ (r^) is likely to be very 
noisy. Thus, differentiation of Ugg^, is not practical. 

We thus consider an alternative approach. Consider a shell with mass dm = Airpr'^dr. Assum¬ 
ing (Uggg) = —2(/)(r), we may write 


dm 


-2vr(uD 


p{r) 


(j){r) 


-dr. 


(18) 


The mass surrounding the halo is 


GM{< r) — GM{< rs) = j {vl^^{x))T{x)dx 

Jrs 


(19) 


where we introduce the filling function J-{r) = —2'kG p{r)r‘^/<p{r). Here, we integrate the noisy 
function rather than differentiating it. 


We now replace, in equation (1£), the three-dimensional distance r with the projected dis¬ 
tance r^, and the three-dimensional escape velocity Uesc with the measurable ui.o.s.- The second 
replacement implies the introduction of the velocity field anisotropy parameter (3. Equation 
becomes 

GM{< r_L) - GM{< rs) = [ {vl^^ {x))J^{x)dx (20) 

Jrs 

where 

3 - 2P{x) 


Tix) =T{x)- 


I- I3{x) 


( 21 ) 


Knowledge of the filling function T allows the estimation of the system mass for any r_i_ G 
(0,oo). However, we need to know T precisely only for r±_ < rs- For rj_ > rs we may consider T 
as a free function, generally slowly varying when r^jrs G (1, 3), where we expect to apply equation 
(1^ . In fact, the assumption of a spherically symmetric p{r) leading to equation (|^) is a very 
crude approximation at best. Thus, the complete expression in equation (^) for T might not 
be a robust representation of the filling function. On the other hand, we will see in §4.2.3 that 
equation (20) still holds for reasonable models of T. In practice, we can estimate the mass within 
rs by applying the virial theorem or by assuming hydrostatic equilibrium of the X-ray emitting gas. 
Outside rs these methods break down, and we can use equation 
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Reasonable assumptions about the behavior of p and (j) at large radii provide a model for T. 
For example, if we neglected the integral in equation (1^), we would have p oc r~^, cp oc r~^, and 
J- = const. A universal density profile suggested by Navarro, Frenk & White (1995) for dark matter 
halos is 


p{r) = 


Pori 


r{r + TsY 

where is some scale length. This profile yields 

~ 1 

nr) = 


(r + TsY 21n(l + r/vs)' 


( 22 ) 


(23) 


This profile is a good fit for r ^ rs only. Gravitational lensing observations show that this model 
reproduces real cluster profiles within 1 MpcQ (Tormen, Bouchet, &: White 1996). However, the 
model may be a poor approximation in the cluster outskirts. 


To obtain we can assume /3(r) ~ const ~ 0.6 — 0.8. Despite this crude approximation to 
jr, we now show that, for suitable choices of the method outlined here can lead to remarkably 
good mass estimates for A^-body systems. 


4. iV-BODY SIMULATIONS 

We now use A-body simulations to illustrate how random motions in the outskirts of dark halos 
invalidate application of the spherical infall model to estimate Hq. We also show that measurement 
of the amplitude of the velocity field does provide a method of obtaining the enclosed mass with 
better than 30% accuracy for massive halos. 


4.1. Simulation Parameters 


As typical hierarchical clustering scenarios, we consider three Cold Dark Matter (CDM) uni¬ 
verses with [Do,Dao] = [1.0,0.0], [0.2,0.0], [0.2,0.8], and h = 0.5. Initial conditions are generated 


with the COSMOS package developed by E. Bertschinger (1995) which perturbs initial particle 
positions and velocities from a grid according to the Zel’dovich approximation. We produce a re¬ 


alization of the Gaussian random field of the initial density perturbations with the Bardeen et al. 


(1986) power spectrum containing the transfer function for adiabatic fluctuations and negligible 
baryon density. We normalize the power spectrum with ag = 1.0 at the present time, where ug is 
the rms matter density fluctuation in spheres of radius 16 Mpc. This value of ug is between the 
values required to fit the observed abundance of local clusters in critical and flat low-density CDM 
universes ( [White, Efstathiou, fc Frenk 1993| ). We simulate a 50^ Mpc^ periodic volume with 64^ 
particles. 


■^All distances in this paper are for h = 0.5. 
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The number density of Abell clusters with richness i? > 1 is ~ 10“® Mpc“^ (e.g. Scaramella et 
al. 1991| ); thus we need constrained initial conditions to obtain a rich cluster within the simulation 
box. We generate constrained random fields with the algorithm of Hoffman & Ribak (1991) in the 
implementation of van de Weygaert & Bertschinger (1996). We are interested only in the presence 
of a rich cluster at the center of the box. Thus, we specify only the height 5p of the local density 
maximum. Specifically, we choose 5p = 3a{RG), where a^{RG) is the variance of the density field 
smoothed with a Gaussian filter of radius Rg = 4 Mpc. 


We integrate the equations of motion with a particle-mesh code ( Hockney &: Eastwood 1981 
Efstathiou et al. 1985| ) using 128^ mesh points, a seven-point finite-difference approximation to 
the Laplacian, a cloud-in-cell density assignment, a leap-frog integrator, and an energy conserving 
scheme to compute the force. CDM universes have an effective spectral index n ~ —1 on cluster 
scales. Thus, we use the scale factor a as the time variable (Efstathiou et al. 1985). We also 
use this variable for the open models, though this choice is actually appropriate only for flat 
universes. Simulations obtained with different integration variables a“ with a £ [0.5,1.5] do not 
yield appreciably different results. Simulations run from a ~ 0.02 to the present time a = 1 
with ~ 700 timesteps. The mesh and the poor momentum conservation of our code produce force 
anisotropies which prevent the Layzer-Irvine cosmic energy equation from being satisfied to better 
than ~ 10% over the entire simulation, even for timesteps smaller than the timestep we used. 
However, the computed force follows the Newtonian value accurately for distances roughly twice 
the mesh cell size, ~ 0.4 Mpc in our simulations. 

Our simulations have two main shortcomings: (1) poor spatial resolution and integration 
accuracy, and (2) small box size. Despite the first problem, we will see in the next subsection that 
for distances larger than the cell size, our halos roughly reproduce the density and velocity fields 
obtained in simulations with larger dynamical range and higher accuracy (Tormen et al. 1996). 
At any rate, we are currently running another set of simulations with the AP^M code (Couchman 
1991) made available by Hugh Couchman and collaborators (Pouchman, Thomas, &: Pearce 1995) 


and we will report on these results in a forthcoming paper. The second problem, namely the small 
box size, implies less tractable errors. 

First, periodic boundary conditions produce artificial tidal forces from replicas of structures 
within the simulation box. The constrained random field clearly aggravates this problem. Artificial 
tidal forces erroneously increase the random motions in the outskirts of halos; they thus artificially 


worsen the prediction of the spherical infall model. Gelb &: Bertschinger (1994b) suggest that box 
sizes L > 50 Mpc are necessary to provide reliable results on clustering. However, in L ~ 50 Mpc 
box simulations, individual halo dynamics is almost unaffected because the artificial tidal forces 
are sufficiently suppressed (pelb Sz Bertschinger 1994^ ). Thus, we expect tidal fields to be a minor 
problem in our models. 


Secondly, we omit contributions from Fourier components of the primordial density perturba¬ 
tion field with wavelengths larger than the box size. Thus, we underestimate both the amplitude 
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of clustering and the amplitude of the peculiar velocities. Corrections to the evolved state of the 
simulation ( Tormen &: Bertschinger 1996| ) or correction during its evolution ( Pole 19^ ) alleviate 
the problem. The missing clustering would probably increase the merger rate and therefore the 
random motion in the outskirts of halos, thus worsening the predictive power of the spherical infall 
model. The velocity field is affected mainly in the bulk velocity. Random motions increase by only 


a few percent when one increases the size of the box from 100 to 800 Mpc (Tormen fc Bertschingi 


mr 


19961) 


In summary, we increase the random motions because of poor spatial resolution and integra¬ 
tion accuracy, but we decrease it because of missing large-scale power. We plan further work to 
investigate how these problems affect our results. 


4.2. Results 

We now describe the results of our simulations. First, we consider the density and velocity 
fields of the most massive halos (§4.2.1). In §4.2.2, we examine the velocity field within the infall 
regions. We compare the simulations with the mass estimation method in §4.2.3. 


4-2.1. Halo Identification 


We identify halos using a generalization of the friends-of-friends algorithm ( Barnes fc Efstathiou 
1987 ) with a linking length equal to 0.1 Mpc ~ 0.13 times the interparticle distance, and a critical 
number of neighbor particles W = 10 = 1 yields the classical friends-of-friends algorithm). 

These parameters assure a halo overdensity 6 ~ 10^ —10^ with respect to the background. However, 
this point is of secondary importance; we are interested only in the determination of the center of 
mass of the halo in order to study the velocity field of the outer regions. Linking lengths between 
0.08 and 0.16 Mpc move the center of mass of the most massive halos less than 1%. The center of 
mass velocity suffers larger oscillations. However, to suppress this problem we redefine the cluster 
velocity from all the particles within the virial radius defined below. 

In the outskirts of halos, halo asphericities have little effect on particle velocities, because 
equipotential surfaces are rounder than equidensity surfaces. For example, van Haarlem et al. 
(1993) find that at distances larger than a few megaparsecs from the halo center, the force differ¬ 
ence between a spherical distribution and a significantly ellipsoidal distribution is ^ 15%. These 
authors also show that asphericities can blur the caustics expected in the infall model. However, 
the disagreement between the spherical infall model and the simulations is more severe than the 
departures induced by asphericity (see van Haarlem &: van de Weygaert 1993D. 


We intend to compare the amplitude of the velocity field with the gravitational potential. Thus, 
we can assume spherically symmetry, even though halos are not spherically symmetric. Around the 
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halo center of mass we define spherical shells at equally spaced logarithmic intervals of distance r 
from the center. Spherical shells containing the same number of particles yield the same results. 


In order to compare our simulations with the spherical infall model, we need to determine 
the region where the model is valid. Consider the halo density at virialization 
where pm is the halo density at maximum expansion, is given by equation ( 0 ), and r] = 

2^7Aoa;max(<5^^)/^o(l + ^i) = A/dvrGpm- We obtain 


Ph _ 2nAo 
p{a) Oo 


(24) 


where p{a) is the mean density of the universe at time a. Note that in Hao 
and equation ( |^) reduces to 

Ph _ 

p{a) Qo ■ 


0 universes 2 flA 0 /p = 


(25) 


In particular, for the Einstein-de Sitter universe the overdensity of the virialized halo scales as 
(I + Zf)^, where Zf is the formation redshift of the halo ( [White 1996| ). With equation (24) we 
can define the radius of the halo rs when the average overdensity 5{r) within the radius r satisfies 
1 + 5 = ph/p{a). Thus, we expect the spherical infall model to be valid when r > rs- 


We now consider the halo density and velocity fields. We consider the most massive halos in 
the simulations at the present time a = 1, when they are far enough in time from major mergers 
which formed the final halos. Strongly unrelaxed states affect the density and velocity fields of the 
halos ( Tormen et al. 1996 ) but have little effect on the predictive power of the escape velocity as 
we show in the next subsection. 


Here we examine the halo formation in the three investigated cosmological scenarios on the 
basis of simulations starting with the same seed for the random number generator. Thus, we obtain 
roughly the same evolutionary history for the three halos. We ran another set of simulations with 
different seeds for the random numbers. These simulations provide similar results. Note that we 
use a fixed number of particles within the simulation box; thus, the halo in the flat universe with 
a final mass ~ 7 ■ Mq is approximately five times more massive than the most massive halos 
in the open universes. 


At a = 1, the largest halos contain 20,000 particles within ~ 2 Mpc, typically. The 
first row of Fig. 2 shows the circular velocity profile Ucirc(^) = [GM{< r)/r]^/^. We compare the 
simulated profiles with the profile 


= dTrCpor; 


I 


In 1 + 


I 


(26) 


r + Tsj 

derived from the Navarro et al. (1995) density profile (eq. [^). Fig. 2 shows the best fits. The 
agreement is acceptable, even though we fit the profile up to r = 5rs instead of r ~ for which 
equation (^) is expected to hold ( Navarro et al. 1995 ). If we fit only the range [0.1, the 
agreement improves slightly, though our poor resolution at small radii is more apparent. 
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The second row of Fig. 2 shows the anisotropy parameter /3(r) = 1 — {v^)/2{v^) (eq. II)- 
When r rs, f3 increases from ~ 0.5 to ~ 0.6 — 0.8, indicating predominantly radial motion. At 
r ~ 5rs, f3 drops to ~ 0.2 — 0.3 because of the presence of another halo. Note that in the flat 
universe, random motions dominate the velocity field more strongly than in the open universes at 
large r. This effect is not the result of the different cosmology, but of the different local dynamics 
and of the larger mass of the halo, which implies a larger accretion rate (see e.g. Lacey fc Cole 
1994 ; Manrique Sz Salvador-Sole 19961) . For comparison. Fig. 2 shows f5{r) for a ~ 2 • IO^^Mq halo 


within the flat universe simulation box (dashed line). The frequency of radial orbits is larger than 
for the more massive halo. 


Our results agree acceptably with those obtained by Tormen et al. (1996) with higher resolution 
simulations of a P{k) oc k~^ flat universe. Note that we reach only a resolution of ~ O.lr^, compared 
with the ~ O.Olr^ resolution of Tormen et al. (1996). 


4-2.2. Velocity Field within the Infall Regions 

We now consider the prediction of the spherical infall model. We project the halo along a 
random direction and consider the amplitude of the velocity field 2Av{r^) = Umax — i^min where 
'f^max and Umin are the maximum and minimum line-of-sight velocity at projected distance r± from 
the center of the halo. Fig. 3 shows Av/Hqt^ in our simulations. This quantity should agree 
with the spherical infall prediction for r±_ > rs (dashed line). We compute this prediction using 
complete three-dimensional information to evaluate the overdensity and then to derive the expected 
velocity field. The model systematically underestimates the amplitude of the velocity field. Clearly, 
estimates of Hq based on spherical infall always overestimate the actual value. 

In contrast, agreement with the escape velocity = —24>{r±) (solid line) is excellent. We 
compute the escape velocity through the discrete version of equation ([l^ . We replace the upper 
limit of integration with a maximum radius rmaxj usually in the range 5 — lOr^. Our results 
are insensitive to this parameter; Umax < usually underestimates the true potential. When 
^'max > lOr^ the spherical assumption breaks down severely because of the presence of other halos. 
We hnally correct for the velocity field anisotropy (eq. [pTf ): 

The agreement is remarkable. It indicates that the gravitational potential, i.e. the local dynamics, 
determines the amplitude of the velocity field in the outskirts of halos. The spherical infall model 
is a poor predictor of the velocity field. The global density of the universe apparently plays no role. 

Fig. 3 also shows that the agreement holds regardless of the dynamical state of the halos. For 
comparison, we show the halos at the present time (lower row), when they are approximately in 
equilibrium, and the halos right after the merging of the two halos of comparable size which formed 
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the final halos (upper row). The escape velocity, namely the gravitational potential, still follows 
the velocity field correctly. 


4-2.3. Estimate of the Halo Surrounding Mass 

We now apply the method outlined in §3 to estimate the mass surrounding the virialized region. 
We suppose that we know the virial radius r^, the particle velocities along the line of sight, and the 
mass within rs, M{< rs). For real systems, we can estimate M(< rs) with X-ray methods or with 
the virial theorem. In future work we plan to investigate how the uncertainty in the virial mass 
affects the estimate of the mass within the infall regions. 

We use equation (^), assuming two different filling functions E, namely Ei = 1/2, and 
E 2 = [ln(l -|- ar_\_/rf)\~^. T\ assumes that the entire gravitational potential originates within r±. 
E 2 is equation with a free parameter a = r^/r* ~ 3 — 10, typically. For E 2 , we have assumed 
r^/(r-|-rs)^ ~ 2(1—/3)/(3 —2/3). This assumption is reasonable when 1 ^ r/r^ ^ 5 and (3 ~ 0.6 —0.7. 
Results obtained with a G (3,10) do not differ appreciably. 

Fig. 4 shows the ratio between the estimated mass within the radius r, Mest(< r) and the 
actual mass M(< r). Equation (p0|) is apparently a good mass estimator, at least up to r ~ brs, 
where we start observing the velocity field of another halo. For r < the difference between 
the estimated and the actual mass is always < 16%, and is ~ 10%, on average. Fig. 4 shows both 
the “unrelaxed” (upper row) and “relaxed” state (lower row; see Fig. 3). For the relaxed state the 
mass estimate improves. However, the difference from the unrelaxed state is not large. 

If we do not know M(< r^) the agreement between the estimated mass and the true mass 
worsens, particularly for unrelaxed states. Fig. 5 shows the halo mass estimated with equation 
(|20t) when we integrate from r_L = 0. It is remarkable that Ei still works reasonably well for r± < rs 
when the clusters are relaxed (lower row), showing that the filling function Ei is more robust than 
E 2 . We expect this result; when r± < rs our hypotheses leading to E 2 break down. When the 
clusters are unrelaxed (upper row) the estimation method introduces large errors. These results 
imply only that both filling functions are inadequate; the velocity field is still well described by the 
escape velocity (Fig. 3). 

Fig. 6 shows how our mass estimator behaves statistically. For each model, we separate the 
most massive from the least massive halos, regardless of their dynamical state. For the flat universe, 
we consider halos with M(< rs) > IO^^Mq and IO^^Mq < M{< rs) < IO^^Mq. For the open 
universes, we consider halos with M{< rs) > 2 • IO^^Mq and IO^^Mq < M{< rs) < 2 ■ IO^^Mq. 
We thus have roughly the same number of particles within each halo for the low or high mass range 
in both the flat and open universes. Fig. 6 shows the median of the halo mass profiles; error bars 
indicate the upper and lower quartile of all the profiles at each position r±frs. 

The method yields better estimates for more massive halos (upper row). The agreement for 
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less massive halos worsens because tidal fields caused by neighboring larger halos severely perturb 
the velocity fields of the smaller halos. However, these small halos contain only a few thousand 
particles within ~ a few times the cell size. Thus, numerical artifacts may also play a significant 
role. 


Fig. 6 shows that our method of estimating the mass of non-virialized regions yields better than 
30% statistical accuracy for halos with virial mass M(< rs) ^ 0.5 — 1.0 ■ IO^^Mq. For comparison, 
consider the accuracy of X-ray estimates of the mass of virialized halos. X-body/hydrodynamics 
simulations show that the assumption of an isothermal gas in hydrostatic equilibrium yields virial 
masses with better than 20% accuracy ( Navarro et al. 1995 ; Schindler 1996| ; [Evrard, Metzler" 
I&: Navarro 1996). We note however that gravitational lensing methods, which do not require 


equilibrium assumptions, may disagree with X-ray observations. For example, Wu fc Fang (1996) 


claim that X-ray observations may actually underestimate the mass of real clusters by a factor of ~ 
2. However, gravitational lensing methods suffer systematic errors due to non-spherical symmetry, 
projection effects, and substructure ( [Miralda-Escude &: Babul 1995 ; Bartelmann 199^ ). On the 
other hand, there are examples where X-ray mass estimates do agree with weak lensing estimates 
to within the error limits (see e.g. Squires et al. 1996a; [Squires et al. 1996b| ). We conclude that 
when applied to X-body systems our mass estimation method has an accuracy comparable with 
other widely used methods. 


Application of our method to real clusters introduces non-trivial challanges. First, we do not 
know how reliably galaxies trace the velocity field in these non-linear regions. Secondly, sampling 
effects may introduce large errors, though these effects can be quantified. Even if we could overcome 
these problems, it is not clear how reliably we can determine the amplitude of the velocity fields in 
cluster infall regions; the caustics may not be very apparent (e.g. [van Haarlem et al. 199^ ). We 
plan to investigate these issues in future work. 


5. CONCLUSION 


In redshift space, galaxies around clusters should appear within regions with a characteristic 
trumpet shape ( Kaiser 1987 ). Regos &: Geller (1989) suggested applying the spherical infall model 
to these caustics to constrain the density of the universe. We show that this method generally 
overestimates the actual density parameter Hq because random motions increase the amplitude of 
the caustics (see also Lilje fc Lahav 19^ [White fc Zaritsky 1992[). 


We define the amplitude of the velocity field as half of the difference between the maximum and 
minimum line-of-sight velocity at projected distance r± from the cluster center. We use W-body 
simulations to show that the escape velocity (eqs. jjl^, and @]) describes this amplitude well, with 
r± ranging over two orders of magnitude, from the central region of the halo to its infall regions, 
where particle orbits are mainly radial. Van Haarlem (1992) first noted that the spherical infall 
method overestimates Hq and pointed out that mergers and substructures within the infall regions 
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are responsible for the disagreement. Here we suggest a unifying explanation. 


We show that our interpretation of the amplitude of the velocity field within halo infall regions 
can be applied to estimate the interior mass of halos up to a few virial radii rs from the halo center, 
where the usual equilibrium assumptions do not hold. This estimation technique works because the 
local dynamics depends more strongly on the mass of the halo than on the global properties of the 
universe. This deduction agrees with the recent suggestion by White (1996) and p^avarro, Frenk 


fc White (1996) that halos have a universal density profile, basically independent of the cosmology, 
with a characteristic density depending on the formation time of the halo (see also our eq. (p^ ) 
which ultimately depends on the halo mass. In the literature, the dependence of the density profile 
on the halo mass has been neglected, with attention focused on the dependence on Hq and on the 
power spectrum P{k). The limited overdensity range examined in these numerical studies explains 
why Wbody simulations appear to show a dependence of the halo density profile on the underlying 
cosmology (e.g. Crone, Evrard, &: Richstone 1994]) . 

If we know the mass M(< r^) within the virial radius our mass estimation method can 
estimate the mass up to several rs with a ^ 30% uncertainty on average, at least for halos with 
mass M(< rs) ^ 0.5 — 1.0 • IO^^Mq. If cluster galaxies are unbiased tracers of the gravitational 
potential, we can estimate the mass within the outskirts of observed systems on the basis of redshift 
data alone. However, we must investigate how sampling affects the result, and we must explore 
how accurately we need to know the mass within the virial radius. Last but not least, we must find 
a reliable method of determining the caustics and extracting the escape velocity function Uesc(^±)- 
This problem is not trivial because the caustics may not be obvious (e.g. van Haarlem et al. 19931) , 
although they do appear in some cases (see e.g. A3266 in [Quintana, Ramfrez, &: Way 1996 , and 
the Coma cluster in the 15R Survey of peller et al. 19^ ). On the theoretical side, we need a 
reasonable assumption for the filling function T that we can calibrate with Wbody simulations. It 
is reassuring that the simplest assumption we investigate here, T = const, works reasonably well. 

We plan to pursue this approach by acquiring dense redshift samples in the infall regions of 
nearby rich galaxy clusters in order to estimate masses on scales ^ 10 Mpc, where linear theory 
breaks down and where galaxy systems are not yet in virial equilibrium. 


We thank Ed Bertschinger for developing the COSMICS package, the cosmological initial con¬ 
dition generator, and for making it available to the scientific community. The COSMICS package 
is funded by the NSF grant AST-9318185. We also thank an anonymous referee for several con¬ 
structive suggestions which improved the presentation of our results. This research is supported in 
part by NASA grant NAGW-201 and by the Smithsonian Institution. 
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Fig. 1.— Amplitude of the velocity field within the infall region of a spherical perturbation 
in an otherwise uniform universe at the present time. The amplitude is in units of the projected 
physical distance r_L from the center of the perturbation; x± is the component perpendicular to 
the line of sight of the local scale factor x. From top to bottom, solid lines are for [no,hlAo] = 
[1.0, 0.0], [0.5, 0.0], [0.1,0.0]. Upper dashed line is for [no,llAo] = [0.5, 0.5]; lower dashed line is for 
[no,nAo] = [0.1,0.9]. 
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[ 1 . 0 , 0 . 0 ] [ 0 . 2 , 0 . 0 ] [ 0 . 2 , 0 . 8 ] 



Fig. 2.— Density and velocity fields of the most massive halo in each A^-body model. Upper row 
shows the circular velocity Vdrc = [GM{< r)/r]^/^. Solid lines are the best fits to equation 
Lower row shows the velocity field anisotropy parameter /3(r) (eq. |l^). The values of [Uoj^Ao] 
are shown over each column. The dashed line in the [Uoj^Ao] = [1-0,0.0] model is for a halo of 
mass M ~ 2 • IO^^Mq, roughly the most massive halo in the open models. 
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Fig. 3.— Amplitude of the velocity field within the infall region of the most massive halo in each 
A^-body model. Solid lines are the escape velocities computed with equation ( [T^ ) corrected for the 
anisotropy parameter with equation (|l^). Dashed lines are the spherical infall predictions which 
hold only for r±_ > rs- Lower row shows the halos at the present time a = 1, when the halos are 
roughly in equilibrium. Upper row shows the halos at earlier times, right after the major mergers 
of the two smaller halos which formed the final halos. 
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Ti/rj r^/rj r^/rj 


Fig. 4.— Ratio of the interior mass Mest{< f) estimated with equation (^) and the true interior 
mass M(< r) of the most massive halo in each A^-body model. We assume that we know and 
the virial mass M(< r^). Bold lines are for the filling function = 1/2. Solid lines are for 
T 2 = [ln(l + Times are as in Fig. 3. 
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[ 1 . 0 , 0 . 0 ] [ 0 . 2 , 0 . 0 ] [ 0 . 2 , 0 . 8 ] 



Fig. 5.— Same as Fig. 4, but we now assume that we do not know the virial mass M(< rs) and 
we integrate equation (^) from = 0. 
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Fig. 6a.— Median mass profiles of halo samples in each A^-body model. Masses are estimated 
with eqnation (p^), assnming we know each M{< rs). Upper rows show the most massive halos: 
M(< rs) > IO^^Mq for the flat model, and M{< rs) > 2- IO^^Mq for the open models. Lower rows 
show the least massive halos: IO^^Mq < M(< rs) < 10^^Mq for the flat model, and IO^^Mq < 
M(< rs) < 2-lO^^M0 for the open models. Numbers of halos in each sample are shown. Error bars 
indicate upper and lower quartiles at each projected distance r±. (a) Filling function = 1/2; 
(b) filling function = [ln(l + 
































































